#include "fermiqcd.h"

void metropolis_multihit(gauge_field &U, 
			 float beta, int n_iter, int n_hits) {
  int iter, parity, hit, mu, ndim=U.lattice().ndim;
  mdp_complex newU; 
  float delta_action;
  mdp_matrix stapleh;
  mdp_site x(U.lattice());
  for(iter=0; iter<n_iter; iter++) {
    for(parity=0; parity<=1; parity++) {
      forallsitesofparity(x,parity) {
	for(mu=0; mu<ndim; mu++) {
	  stapleh=staple_H(U,x,mu);
	  for(hit=0; hit<n_hits; hit++) {
	    newU=exp(2.0*Pi*I*U.lattice().random(x).plain());
	    // CHECK FOLLOWING LINE FOR SIGN AND COEFFICINTS
	    delta_action=beta*real(trace((newU-U(x,mu))*stapleh)); 
	    if(delta_action<U.lattice().random(x).plain()) U(x,mu)=newU;
	  }
	}
      }
      U.update(parity);
    }
  }
}

int main(int argc, char** argv) {
  mdp.open_wormholes(argc,argv);
  int L[]={4,4,4,4}; // the lattice size
  mdp_lattice spacetime(4,L); // in 4D
  gauge_field U(spacetime,1); // 1 here means SU(1)
  // set cold
  set_hot(U);
  for(int i=0; i<10; i++) {
    metropolis_multihit(U,4.0,10,2); // beta=4.0, 10 steps, 2 hits
    mdp << "P=" << average_plaquette(U) << endl;
    // GaugeFixing::fix(U,GaugeFixing::Landau,100,1e-5); 
    // 100 max steps and 1e-5 target precision
    mdp << "P=" << average_plaquette(U) << endl;
    // do whatever you want with it
  }
  mdp.close_wormholes();
  return 0;
}
